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Gravitational  and  magnetic  anomaly  inversion 
using  a  tree-based  geometry  representation 


Raymond  A.  Wildman1  and  George  A.  Gazonas1 


ABSTRACT 

Gravitational  and  magnetic  anomaly  inversion  of  homoge¬ 
neous  2D  and  3D  structures  is  treated  using  a  geometric  pa¬ 
rameterization  that  can  represent  multiple,  arbitrary  poly¬ 
gons  or  polyhedra  and  a  local-optimization  scheme  based  on 
a  hill-climbing  method.  This  geometry  representation  uses  a 
tree  data  structure,  which  defines  a  set  of  Boolean  operations 
performed  on  convex  polygons.  A  variable-length  list  of 
points,  whose  convex  hull  defines  a  convex  polygon  operand, 
resides  in  each  leaf  node  of  the  tree.  The  overall  optimization 
algorithm  proceeds  in  two  steps.  Step  one  optimizes  geomet¬ 
ric  transformations  performed  on  different  convex  polygons. 
This  step  provides  approximate  size  and  location  data.  The 
second  step  optimizes  the  points  located  on  all  convex  hulls 
simultaneously,  giving  a  more  accurate  representation  of  the 
geometry.  Though  not  an  inherent  restriction,  only  the  geom¬ 
etry  is  optimized,  not  including  material  values  such  as  densi¬ 
ty  or  susceptibility.  Results  based  on  synthetic  and  measured 
data  show  that  the  method  accurately  reconstructs  various 
structures  from  gravity  and  magnetic  anomaly  data.  In  addi¬ 
tion  to  purely  homogeneous  structures,  a  parabolic  density 
distribution  is  inverted  for  2D  inversion. 


INTRODUCTION 

Gravitational  and  magnetic  anomaly  inversion  has  applications  in 
many  fields,  including  geophysical  prospecting  and  archeology 
(Oldenburg,  1974).  Fortunately,  exact  forward  models  for  comput¬ 
ing  gravitational  and  magnetic  fields  caused  by  either  infinite  polyg¬ 
onal  cylinders  (2D)  or  arbitrary  polyhedra  (3D)  exist  so  that  compu¬ 
tational  modeling  and  inversion  of  such  structures  is  very  efficient 
(Hubbert  1948;  Won  and  Bevis,  1987).  Several  methods  for  model¬ 
ing  and  inversion  have  been  explored,  such  as  those  based  on  a  data¬ 
base  of  known  masses  (Bullard  and  Cooper,  1948),  linear  splines 


(Murthy  and  Rao  1993),  neural  networks  (Osman  et  al.,  2006,  Os¬ 
man  et  al.,  2007)  continuous  curves  (Abdelrahman  and  Essa,  2005; 
Essa,  2007),  linearization  of  the  nonlinear  integral  equation  (Gao  et 
al.,  2007),  singular  value  decomposition  (Lines  and  Treitel,  1984), 
Fourier  transform  (Mareschal,  1985),  simulated  annealing  (Mundim 
et  al.,  1998),  2D  binary  grid  methods  (Krahenbuhl  and  Li,  2006),  and 
3D  prism  methods  (Rao  et  al.,  1999). 

The  method  presented  here  focuses  on  the  data  structure  used  to 
represent  the  geometry  of  the  inverted  structure.  First  introduced  in 
Wildman  and  Weile  (2007),  the  data  structure  is  a  binary  tree  that  de¬ 
fines  a  set  of  Boolean  operations  performed  on  convex  polygons. 
Each  convex  polygon  is  defined  as  the  convex  hull  of  a  list  of  points. 
This  representation  can  then  use  line  segments  to  approximate  any 
geometry  or  topology.  The  method  can  represent  exactly  any  poly¬ 
gon  and  also  approximate  curved  structures  using  an  arbitrary  num¬ 
ber  of  line  segments.  Also,  multiple  shapes  are  easily  represented  us¬ 
ing  the  tree  structure;  because  the  size  of  the  tree  is  unrestricted,  any 
number  of  shapes  can  be  imaged. 

This  method  has  a  number  of  advantages.  First,  contrary  to  spline- 
based  methods,  the  number  of  points  used  in  any  single  convex  hull 
is  not  restricted.  In  spline  approaches,  using  too  few  points  can  lead 
to  an  inaccurate  geometry.  If  too  many  points  are  chosen,  the  optimi¬ 
zation  can  be  inefficient.  Second,  the  number  of  separate  shapes  is 
unrestricted.  Again,  many  spline-based  methods  require  a  priori 
knowledge  of  the  number  of  separate  geometries,  and  the  possible 
intersection  between  shapes  can  be  difficult  to  handle.  Third,  the  size 
of  the  search  space  scales  up  only  as  the  geometry  becomes  more 
complicated.  In  contrast,  pixel-based  methods  rely  on  a  discretiza¬ 
tion  of  the  region  that  can  be  too  coarse  or  too  fine.  Fine  grids  are  ca¬ 
pable  of  a  higher-resolution  image;  however,  the  search  space  is 
much  larger  at  the  outset.  Fourth,  Fourier-based  methods  cannot  re¬ 
produce  nonsmooth  structures.  The  use  of  line  segments  allows  for 
accurate  representation  of  comers,  and  because  their  length  can  be 
arbitrarily  small,  curves  are  also  well  approximated.  Finally,  the 
method  easily  extends  to  three  dimensions.  Every  geometric  concept 
used  has  an  obvious  3D  analog.  (See  3D  examples  in  the  results  sec¬ 
tion;  however,  the  algorithm  is  not  described  in  detail  because  it  is 
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essentially  identical  to  the  2D  algorithm.)  In  total,  our  method  makes 
no  assumptions  regarding  the  number  or  complexity  of  the  objects  to 
be  imaged. 

Previously,  a  genetic  algorithm  (GA)  was  used  with  the  above  ge¬ 
ometry  representation  for  the  electromagnetic  inversion  of  conduct¬ 
ing  cylinders  (Wildman  and  Weile,  2007).  Although  the  GA  ap¬ 
proach  showed  the  viability  of  the  method,  it  required  a  significant 
number  of  function  evaluations;  subsequently,  a  greedy  search 
method  and  a  combined  GA/local  search  method  were  implemented 
to  improve  performance  (Wildman  and  Weile,  2008).  The  method 
discussed  here  forgoes  any  stochastic  optimization  and  adapts  a 
Newton-like  local  search  method,  giving  a  simplified  and  more  effi¬ 
cient  approach.  Ultimately,  the  method  uses  local  search  combined 
with  geometric  operations  discussed  in  previous  work. 

The  local  search  is  used  to  optimize  two  quantities:  affine  trans¬ 
formations  that  are  applied  separately  to  each  convex  polygon  in  a 
tree  and  the  individual  points  on  the  convex  hulls  of  each  polygon. 
Each  optimization  scheme  has  advantages.  The  affine  transforma¬ 
tion  is  lower  dimensional,  leading  to  more  efficient  convergence  of 
the  local  search.  Approximate  size  and  location  of  geometry  can  be 
determined  using  this  optimization.  Optimizing  the  individual 
points  is  more  efficient  for  recovering  detailed  shape  information. 
The  overall  optimization  algorithm  proceeds  in  two  stages,  only  us¬ 
ing  higher  dimensional  searches  after  other  options  have  been  ex¬ 
hausted.  The  result  is  an  efficient  search  algorithm  that  typically  re¬ 
turns  results  with  low  relative  errors  in  the  gravitational  or  magnetic 
field  (on  the  order  of  10~3  to  10-4  as  measured  in  the  €2-norm  for 
simulated  targets)  using  on  the  order  of,  depending  on  the  complexi¬ 
ty  of  the  target,  103  to  104  function  evaluations. 

In  the  next  section,  the  geometry  optimization  method  is  detailed. 
The  geometric  data  structure  is  reviewed,  and  a  local  optimization- 
based  method  for  inversion  is  discussed.  The  numerical  results  sec¬ 
tion  explores  the  use  of  gravitational  and  magnetic  data  as  well  as 
measured  data  from  the  literature.  The  final  section  discusses  con¬ 
clusions  and  future  work. 

GEOMETRY  OPTIMIZATION 

This  section  details  the  geometry  encoding  and  optimization 
scheme.  In  the  first  subsection,  the  key  data  structure  is  reviewed  and 


Figure  1 .  Convex  hull  of  a  set  of  points. 


discussed.  The  second  subsection  details  the  optimization  algorithm 
which  is  a  combination  of  a  Newton-like  search  method  and  several 
geometric  operations.  The  objective  function,  which  relates  the  error 
in  a  potential  solution  to  some  synthetic  or  measured  data,  is  present¬ 
ed  in  the  third  subsection.  Next,  the  gradient-based  search  method  is 
reviewed  and  discussed  in  context  with  the  given  objective  function 
and  geometry  representation.  Finally,  the  algorithm  is  reviewed  in 
total  through  an  example. 

Geometry  representation 

The  basic  principle  behind  this  geometry  parameterization  is  the 
use  of  convex  polygons  to  build  more  complex  geometries  and  to¬ 
pologies.  Though  unnecessary,  convex  shapes  are  ideal  because  they 
can  be  represented  easily  as  the  convex  hull  of  a  list  of  points.  The 
convex  hull  of  a  set  of  points  is  defined  as  the  smallest  convex  set  that 
contains  all  points  or,  equivalently,  the  intersection  of  all  half-planes 
that  contain  the  points  (de  Berg  et  al.,  2000).  Figure  1  gives  an  exam¬ 
ple  of  the  convex  hull  (solid  line)  of  a  set  of  points  (shown  as  black 
dots). 

Any  list  of  three  or  more  noncollinear  points  generates  a  valid 
convex  hull,  making  self-intersections  impossible.  A  mathematical 
definition  of  the  convex  hull  (CH)  of  a  set  of  N  points,  X 
=  {pu) ,  ...  ,  p(w)},  is  given  by 

N  N 

2  «/P<')|«;  &  0  ,  V  i ,  2  «,  =  1  |  ■  (1) 
i  =  1  i  =  1  J 

Though  this  strict  definition  defines  the  set  of  (infinitely  many) 
points  contained  in  a  convex  hull,  typical  algorithms  return  the  sub¬ 
set  of  points  that  are  located  on  the  vertices  of  a  convex  hull  (as 
shown  in  Figure  1  and  denoted  as  matrix  P  throughout)  rather  than 
the  set  of  coefficients  a,-. 

In  this  work,  line  segments  are  used  to  connect  points,  giving  a  lin¬ 
ear,  polygon-based  approximation  of  target  geometries.  This  is  not  a 
limitation  of  the  method  because  splines  or  Bezier  curves  could  also 
be  used  (Farin  and  Hansford  2000;  Mortenson,  1999).  Additionally, 
the  number  of  points  on  a  convex  hull  is  not  restricted  in  any  way,  so 
arbitrarily  small  line  segments  can  be  used  to  approximate  curved 
structures. 

More  complex  structures,  such  as  multiple  and  concave  shapes, 
can  be  generated  by  combining  convex  shapes.  One  way  of  combin¬ 
ing  convex  polygons  is  to  generate  expressions  involving  Boolean 
set  operations  such  as  union,  subtraction,  and  intersection  applied  to 
convex  polygon  operands.  Moreover,  mathematical  expressions  can 
be  represented  as  binary  trees,  which  can  be  manipulated  with  an  op¬ 
timization  scheme.  Here,  the  internal  nodes  of  the  tree  represent 
Boolean  operations  and  the  leaf  nodes  represent  convex  polygons. 

One  difference  between  this  work  and  Wildman  and  Weile  (2007) 
is  that  only  union  operations  are  used.  Figure  2  gives  an  example  of 
this  scheme.  Consider  a  Boolean  expression  (shown  in  Figure  2a  as  a 
tree)  of  the  form 

CH(X)  UCH(A)  UCH(O),  (2) 

where  X ,  A,  and  0  represent  three  separate  sets  of  points  (shown  in 
Figure  2b,  with  dashed  lines  representing  the  convex  hulls)  and  U 
represents  union. 

The  first  step  in  generating  the  geometry  defined  by  expression  2 
is  to  compute  the  union  of  the  convex  polygons  represented  by 
CH(  X )  and  CH(  A),  as  shown  in  Figure  2c.  In  this  case,  the  polygons 
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are  separate,  so  computing  the  union  is  trivial.  Next,  the  union  of  the 
result  of  the  previous  step  (shown  in  Figure  2d  as  the  shapes  with  Os 
on  the  vertices)  and  the  convex  polygon  represented  by  CH(  0 )  is 
computed.  In  Figure  2d,  the  result  of  the  union  operation  is  shown  as 
the  solid  line,  and  the  operands  are  shown  as  dotted  lines.  Figure  2e 
illustrates  the  result  of  expression  2. 

Evaluating  geometries  such  as  expression  2  requires  several  com¬ 
putational  geometry  algorithms.  First,  there  are  a  number  of  convex 
hull  algorithms;  should  structures  with  a  large  number  of  points  be 
necessary,  0(n  log(n))  (in  the  total  number  of  points  n )  algorithms 
exist.  Computing  Boolean  operations  requires  an  algorithm  known 
as  map  overlay  (de  Berg  et  al.,  2000).  This  algorithm  computes  a 
partitioning  of  the  plane  based  on  the  input  polygons.  Intersection 
points  and  separate  regions  are  identified  so  that  Boolean  opera¬ 
tions  can  be  performed.  This  algorithm  also  can  be  computed  in 
0(n  log(n)  +  I  log(n))  time,  where  the  computation  time  depends 
on  the  number  of  intersections  I  (de  Berg  et  al.,  2000).  In  this  imple¬ 
mentation,  the  Computational  Geometry  Algorithms  Library  was 
used  for  all  3D  geometry  operations. 

Ultimately,  this  scheme  applies  an  optimization  process  to  the 
data  structure  defined  above.  The  data  structure  can  be  summarized 
as  a  set  of  lists  or  arrays  containing  2D  (or  3D)  points  stored  in  the 
leaf  nodes  of  a  tree.  The  internal  nodes  of  the  tree  define  some  Bool¬ 
ean  combination  of  the  convex  polygons  that  each  point  list  repre¬ 
sents.  An  optimization  scheme  applied  to  this  data  structure  must 
then  operate  on  individual  points  and  the  actual  structure  of  the  tree. 
A  numerical  optimization  scheme  is  applied  to  the  convex  shapes 
(directly  optimizing  the  coordinates  and  geometric  transforma¬ 
tions),  and  decisions  in  the  shape  of  the  tree  are  made  based  on  the 
performance  of  several  runs  of  the  optimizer. 

Optimization  algorithm 

Given  the  possible  complexity  of  the  geometry  representation,  an 
optimization  scheme  must  be  designed  to  construct  a  target  geome¬ 
try  efficiently.  To  that  end,  the  local  optimization  method  used  here 
gently  increases  the  dimension  of  the  search  space  by  optimizing  dif¬ 
ferent  geometric  transformations  and  carefully  controlling  the  size 
of  the  tree.  The  method  proceeds  in  two  stages:  a  split  stage  deter¬ 
mines  the  approximate  size,  location,  and  geometric  complexity  of 
the  target  and  an  optimize  stage  gives  a  more  accurate  image.  Each 
iteration  of  the  split  stage  divides  the  leaf  nodes  of  the  tree  into  the 
union  of  two  separate  convex  polygons  along  one  or  more  angles. 
The  best  choice  is  saved  and  optimized.  Once  the  change  in  the  ob¬ 
jective  function  is  sufficiently  small,  the  optimize  stage  applies  high¬ 
er-dimension  searches  to  the  geometry  to  model  the  target  more  ac¬ 
curately. 

The  optimization  algorithm  used  here  is  the  Broyden-Fletcher- 
Goldfarb-Shanno  (BFGS)  local  search  method  (Dennis  and  Schna¬ 
bel,  1996).  BFGS  is  a  hill-climbing  method  similar  to  Newton’s 
method.  It  differs  from  Newton’s  method  in  that  it  uses  a  specific 
rank-two  update  to  the  Hessian  matrix  of  approximated  second  de¬ 
rivatives.  Because  the  search  space  is  multidimensional,  line  search 
is  used  at  a  fixed  number  of  search  directions  to  find  successive  mini¬ 
mized  vectors.  Throughout  this  paper,  the  term  iteration  refers  to  a 
single  loop  through  a  stage  of  the  global  scheme,  not  separate  search 
directions  within  BFGS. 


The  dimension  of  the  space  searched  by  BFGS  can  be  controlled 
by  optimizing  various  affine  transformations,  including  scaling  and 
translation.  An  affine  transformation,  essentially  a  linear  transfor¬ 
mation  and  a  translation,  can  be  expressed  as 

x*  =  Ax  +  b,  (3) 

where  A  is  a  square  matrix  representing  the  linear  transformation 
and  b  is  a  column  vector  representing  the  translation.  When  applied 
to  a  set  of  2D  points,  expression  3  can  be  rewritten  as 

P*  =  A(P  -  C)  +  B  +  C, 

B  =  b[ll ...  l]lx„> 

C  =  c[ll ...  l]lx„,  (4) 

where  A  is  a  2  X  2  matrix  representing  a  linear  transformation,  b  is  a 
2D  column  vector  representing  a  translation  vector,  c  is  a  2D  column 
vector  representing  the  center  of  mass  of  the  set  of  points,  and  P  is  a 
2  X  n  matrix  representing  the  components  of  the  n  points.  (For  three 
dimensions,  the  dimensions  of  all  matrices  and  vectors  are  increased 
to  3  X  3  and  3  X  1,  respectively.) 


Figure  2.  Decoding  process,  (a)  Tree  structure,  (b)  Point  lists  with 
their  convex  hulls,  (c)  Union  of  two  separate  convex  hulls,  (d)  Union 
with  remaining  shape,  (e)  Final  result. 
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Geometric  transforms  can  be  applied  (and  optimized)  by  selecting 
specific  parts  of  an  affine  transformation  given  by 


A( $XX  '  $XZ  >  $zx  I  '''zz) 


$xx  sxz 

-  $zx  ^  Z,Z  - 


b (tx  ,  tz) 


Scaling  in  two  dimensions  involves  optimizing  the  quantities  sxx  and 
s~,  with 


^i^XX  ’  ‘^Zz) 


o 


(6) 


Similarly,  translation  in  two  dimensions  involves  optimizing  the  ele¬ 
ments  of  b,  tx ,  and  tz. 


1  0 

t„ 

_0  1_ 
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Combining  scaling  and  translation  gives  a  four-parameter  optimiza¬ 
tion  problem: 
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Jz  - 

P  = 


(10) 


Again,  the  initial  guess  used  with  BFGS  is  the  original  set  of  points 
with  some  small  random  variable  added. 

Figure  4  shows  a  flow  chart  of  the  overall  optimization  scheme. 
The  method  is  broken  into  three  stages:  initialize,  split,  and  opti¬ 
mize.  Each  stage  consists  of  a  loop  over  a  set  of  optimizations,  each 
with  a  certain  termination  criterion.  Termination  criteria  can  be  a 
fixed  number  of  iterations,  achievement  of  a  fixed  absolute  error,  or 
stagnation  measured  as  an  insufficient  change  in  the  error. 


Initialize 

The  initial  geometry  is  always  a  centered  rectangle  with  x  and  z  di¬ 
mensions  equal  to  half  the  specified  bounds  of  the  region.  This  step  is 
denoted  as  the  Rectangle  block  in  Figure  4.  Next,  the  initial  geome¬ 
try  is  optimized  using  the  scaling  and  translation  transforms  given  in 
equations  6  and  7.  The  two  optimizations  are  repeated  a  set  number 
of  times;  typically  two  or  three  iterations  are  sufficient.  Using  only 
scaling  and  translation  optimizations  gives  approximate  size  and  lo¬ 
cation  information. 


Split 


Figure  3  shows  an  example  of  an  affine  transformation  applied  to 
a  rectangle,  characterized  by 


1.4 

-  0.1 " 

5 

A  = 

.o.i 

1.2 

,  b  = 

_  1.844  _ 

where  A  is  unitless  and  b  is  measured  in  kilometers. 

In  addition,  a  separate  transformation  is  applied  to  each  point  list 
in  a  tree.  Given  k  leaf  nodes  in  a  tree,  optimizing  a  scaling  or  transla¬ 
tion  transformation  requires  2k  parameters,  each  combined  scaling 
and  translation  requires  4k  parameters,  and  each  affine  transforma¬ 
tion  requires  6k  parameters.  The  initial  guess  for  each  leaf  node  is 
chosen  as  the  corresponding  identity  transformation  with  some 
Gaussian  random  variable  (with  a  small  variance)  added.  (In  three 
dimensions,  the  corresponding  optimizations  require  3k,  6k,  and  12k 
parameters,  respectively.) 

Individual  points  are  also  optimized.  Given  a  set  of  N  points, 
BFGS  must  optimize  a  set  of  2 N  parameters  (or  3 N  in  three  dimen¬ 
sions)  comprised  of  the  components  of  each  point.  The  set  of  optimi¬ 
zation  parameters  is  given  by 


The  split  stage  determines  the  geometric  and  topological  com¬ 
plexity  of  the  structure.  The  splitting  operation  splits  a  chosen  leaf 
node  along  a  given  line  into  the  union  of  two  separate  nodes  while 
adding  the  intersection  points  of  the  line  and  the  hull  to  the  newly 
created  leaf  nodes.  Including  intersection  points  ensures  that  the  re¬ 
sulting  geometry  is  identical  to  the  original.  Figure  5  shows  an  exam¬ 
ple  of  the  splitting  process:  Figure  5a  shows  the  initial  and  resulting 
tree  structure,  and  Figure  5b  and  c  shows  the  corresponding  geome¬ 
tries.  (The  centerline  in  Figure  5c  is  for  illustrative  purposes;  the 
added  intersection  points  ensure  that  the  resulting  geometry  is  iden¬ 
tical  to  the  initial.) 

The  “split  all”  block  attempts  to  determine  the  best  choice  of  split 
between  all  current  convex  shapes.  Each  convex  shape  is  first  split 
along  0°  and  subsequently  90°  and  is  optimized  using  the  combined 
translation/scaling  given  in  equation  8.  (In  three  dimensions,  the  ob¬ 
jects  are  split  along  the  three  planes  that  compose  the  coordinate 
axes.)  After  all  combinations  of  splitting  are  completed,  the  best  re¬ 
sult  is  saved.  Next,  optimization  of  a  full  affine  transformation  is  ap¬ 
plied  to  the  newly  split  structure;  however,  the  result  of  this  optimi¬ 
zation  is  kept  only  if  one  of  the  termination  criteria  of  the  split  stage 
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Figure  3 .  Affine  transformation  applied  to  a  rectangle,  (a)  Initial  shape,  (b)  Affine  transformation  applied. 
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is  achieved.  The  termination  criteria  used  here  are  the  reduction  be¬ 
low  an  error  of  1 0  " 2  or  the  difference  between  the  initial  error  and  er¬ 
ror  after  splitting  below  1CU2.  The  former  criterion  indicates  that  the 
structure  is  likely  accurate  enough  and  further  subdivision  will  not 
aid  the  final  step  because  the  latter  indicates  that  the  splitting  process 
has  stagnated. 

Optimize 

The  final  stage  simply  optimizes  the  divided  structure  alternating 
between  affine  transformations  (equation  5)  and  points  (equation 
10).  In  addition,  new  points  may  be  added  in  between  the  two  opti¬ 
mization  runs.  In  this  implementation,  a  point  is  added  in  turn  to  the 
midpoint  of  each  line  segment  (or  polygon  centroid  in  three  dimen¬ 
sions)  and  slightly  translated  outward  along  the  normal  direction  of 
the  line  segment.  The  objective  function  of  this  structure  is  evaluat¬ 
ed,  and  the  new  structure  with  the  best  objective  function  value  is 
saved.  Finally,  the  termination  criteria  are  similar  to  the  split  stage, 
though  both  the  error  and  difference  in  the  error  must  be  below  1 0  " 3. 
In  addition,  this  step  can  be  terminated  after  four  or  five  iterations. 


station  location  o„  and  index  arithmetic  is  performed  cyclically.  Sev¬ 
eral  special  cases  exist  (e.g.,  rk  —  0),  details  of  which  can  be  found  in 
Won  and  Bevis  (1987).  Also,  magnetic  anomalies  m,  are  computed 
in  a  similar  fashion.  Finally,  gravity  anomalies  of  polygons  with  par- 


Figure  4.  Optimization  algorithm  flowchart. 


Objective  function 

An  objective  function  must  be  defined  to  measure  the  perfor¬ 
mance  of  the  algorithm  relative  to  some  input  data.  The  objective 
function  is  then  the  relative  error  of  a  potential  solution  compared 
with  simulated  or  actual  data  produced  by  a  target  geometry.  Gravi¬ 
tational  and  magnetic  anomalies  can  be  computed  efficiently  for  po¬ 
lygonal  geometries  such  as  those  produced  by  the  proposed  geome¬ 
try  representation  in  Figure  2  (Won  and  Bevis,  1987). 

For  homogeneous  polygonal  structures  (represented  as  K  line  seg¬ 
ments),  the  vertical  component  of  the  gravity  anomaly  measured  at  a 
station  o,  can  be  expressed  as 

K 

Si  =  2Gp^ZM,  (11) 

k=  1 


where  G  is  the  gravitational  constant,  p  is  the  density  contrast  of  the 
structure,  and  Zk  is  a  line  integral  over  the  Tth  segment  of  the  struc¬ 
ture.  Analytical  expressions  for  the  line  integrals  Zk  are  available, 
and  one  efficient  representation  is  given  by  Won  and  Bevis  (1987): 


Z*(o ,-)  =  A 


{0k~ 


6k+ 1)  +  B  In' 


(12) 


where 


A  = 


(P 


<*+D 


Pri)Pf) 


t  (k  + 1)  (k)\2  ,  /  (*+1)  (k)\ 2 

(Px  -  Px  )  +  (pz  Pz  ) 
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B  = 


Jk+  1) 
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(14) 


A  =  (pf])2  +  ( 

(15) 

/  _(*)  \ 

^  =  tan_1(±j- 

(16) 

The  terms  p[k)  and  p[k>  represent  the  coordinates  of  the  kth  vertex  of  a 
structure  (equation  1 0,  or,  e.g.,  the  result  of  equation  2)  relative  to  the 


Figure  5 .  The  splitting  process,  (a)  Tree  structure,  (b)  Initial  point  list 
with  convex  hull,  (c)  Union  of  two  convex  hulls  at  90°. 
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abolic  density  distributions  and  3D  polyhedra  can  be  computed  as 
given  in  Rao  et  al.  (1994)  and  Singh  and  Guptasarma  (2001). 

For  2D  geometries,  a  set  of  Ns  measurements  is  taken  at  linearly 
spaced  stations  along  the  z  =  0,  0  :£  x  £  ymax  line  segment.  (For  3D 
cases,  a  uniform  grid  of  observation  stations  is  used  to  measure  the 
target  data.)  Although  the  altitudes  of  the  stations  are  assumed  to  be 
constant,  local  topography  can  be  modeled  by  varying  the  station 
height  z.  The  objective  function  is  defined  as  the  € 2  errors  relative  to 
some  measured  or  synthetic  data  for  gravitational  data 


/g  = 


and  magnetic  data 


/m  = 


?•  ~  8 

■’oAL  I 

2.-=ik 


meas|2 


meas|2 


2; 


meas|2 


i  =  ir 


I  n 

2,=iK 


(17) 


(18) 


where  g,  and  g”1'”5  represent  the  gravitational  anomaly  at  the  ith  sta¬ 
tion  of  the  search  agent  (as  computed  by  equation  11)  and  target  ge¬ 
ometry,  respectively,  and  m,  and  m ”eas  represent  the  total  magnetic 
anomaly  at  the  ith  station.  Because  the  magnetic  field  m  is  a  vector, 
some  combination  of  the  magnetic  field  components  m„  my,  and  mz 
must  be  used  to  obtain  a  scalar  objective  function.  The  total  magnetic 
field  is  used  and  is  given  by  (for  2D  cases) 

m  =  sin  {(f>mc)mz  +  sin(<^stl)cos(<^lnc)mt,  (19) 

where  (f>mc  is  the  inclination  angle  of  the  ambient  magnetic  field  mea¬ 
sured  in  degrees  below  horizontal  and  is  the  strike  angle  of  the 
structure  relative  to  magnetic  north  measured  in  degrees.  In  three  di¬ 
mensions,  the  total  magnetic  field  is 

m  =  \/m2  +  m2  +  m2.  (20) 

Finally,  noise  can  be  added  to  synthetic  data  to  simulate  a  physical 
system.  Here,  noise  is  defined  as  an  additive  Gaussian  random  vari¬ 
able  with  zero  mean  and  a  standard  deviation  given  by  some  defined 
signal  to  noise  ratio  (S/N),  measured  in  decibels  (dB).  The  measured 
data  is  then  given  by 

8i  =  8i  +  ,  (T)rms{g  ), 


/z  =  0,  o- =  10's/N/2°,  (21) 

where  g)‘m  is  the  simulated  gravitational  field  at  the  ith  station, 
,  a)  is  a  Gaussian  random  variable  with  mean  /x  and  variance  cr, 
and  rms(-)  indicates  the  root-mean-square  average  of  the  simulated 
data.  A  similar  expression  is  defined  for  the  magnetic  data  m™™. 


BFGS  applied  to  geometry  optimization 

BFGS  forms  the  basis  of  the  proposed  inversion  scheme.  It  is  used 
to  optimize  different  sets  of  geometric  transformations  and  also 
point  locations  directly;  however,  the  basic  search  method  is  the 
same  regardless  of  the  chosen  optimization  parameters.  Again,  most 
steps  shown  in  Figure  4  involve  applying  BFGS  to  the  current  geom¬ 
etry. 

To  reiterate,  BFGS  is  a  Newton-like  local  search  method  that  uses 
an  approximation  of  the  Hessian  matrix  of  second-order  partial  de¬ 
rivatives  to  avoid  costly  (and  possibly  inaccurate)  second-derivative 
computations  (Dennis  and  Schnabel,  1996).  The  method  iterates 
through  a  number  of  search  directions  that  are  determined  by  solving 


H,s*  =  -  V  f(xk)  (22) 

for  the  kth  search  direction  sh  where  Hj.  is  the  approximated  Hessian 
matrix,  xk  is  a  vector  of  the  independent  variables  used  to  compute 
either  fg  or/m,  and  the  gradient  V/ is  computed  via  finite  difference. 

The  contents  of  the  vector  of  independent  variables  depend  on 
which  type  of  optimization  is  being  used  because,  as  discussed 
above  (equations  5-10),  scaling,  translation,  affine  transformations, 
or  point  locations  can  be  optimized.  An  affine  transformation  ap¬ 
plied  to  a  single  leaf  node,  for  example,  is  represented  by 

x  [‘At  5  Az  5  A* 5  A?  ’  fv  5  •  (23) 

In  total,  computing /(x)  involves  applying  the  given  transforma¬ 
tion  to  the  current  tree-based  geometry  (equations  5-8),  evaluating 
the  geometry  (e.g.,  as  in  expression  2),  computing  the  gravitational 
or  magnetic  anomaly  (equations  12-16),  and  computing  the  error 
relative  to  the  target  as  given  in  equations  17  and  18. 

With  a  new  search  direction  in  hand,  the  next  point,  xt+1,  can  be 
found  using  a  line  search.  The  line  search  proceeds  by  finding  an  ak 
that  minimizes  the  linear  model 


/*(«*)  =  f(xk  +  ak  sk),  (24) 

giving  a  new  point 

**+t  =  *k  +  «*S*.  (25) 

Given  a  new  vector,  the  Hessian  matrix  must  be  updated  to  find  the 
next  search  direction.  First,  the  gradient  of  the  function  at  the  new 
point  is  computed  and  the  difference  of  the  gradients  at  the  new  and 
previous  points  is  defined  as 

g *  =  V/(x*+1)  -  V/(x*).  (26) 


The  Hessian  matrix  is  then  updated  using  a  rank-two  approximation, 
given  by 


H 


it+i 


H,  + 


gkgk  H*s*(H*s*)J 


T 

g  *s* 


k 


(27) 


which  is  rank  two  because  it  is  the  combination  of  two  rank-one  ma¬ 
trices  (i.e.,  an  outer  product  of  two  vectors),  each  with  a  different  ba¬ 
sis  vector.  Finally,  the  inverse  of  the  Hessian  matrix  is  required  to 
solve  equation  22  and  can  be  expressed  as 


H 


k+  1 


-  Ht.  +  gkgk 


^k^k  g  k  gk  H/;  g k 

(sjgt)2 


_  H;  lgkSTk  +  skgTkHk  1 
sjg  k 

For  the  initial  approximate  Hessian  matrix  H0,  a  diagonal  matrix 
is  formed  using  values  corresponding  to  the  inverse  squared  of  the 
expected  magnitude  of  the  variables  to  be  optimized.  Expected  mag¬ 
nitudes  for  A  are  chosen  as  one,  and  magnitudes  for  b  can  be  chosen 
as  some  fraction  of  the  total  region  size,  given  by  xmax  and  zmax.  For  an 
affine  transformation,  the  initial  Hessian  matrix  is 


Hq  =  diag[l,l,l,l(cxmax)  2[-](czmax)  2],  (29) 

where  diag[-]  indicates  a  diagonal  matrix  with  the  argument  placed 
along  the  diagonal,  and  c  is  a  constant  typically  chosen  as  0.25.  The 
initial  guess  x0  is  chosen  as  the  identity  operation  for  the  chosen  geo¬ 
metric  transformation  or,  for  direct  point  optimization,  the  set  of 
points  of  the  current  geometry.  In  addition,  a  small  random  value 
chosen  from  a  zero  mean  Gaussian  distribution  is  added  to  the  initial 
guess.  The  variance  of  the  distributions  is  determined  by  the  typical 
magnitudes  discussed  above. 
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For  an  affine  transformation,  matrix  A  is  chosen  as  the  identity 
matrix  and  the  elements  of  vector  b  are  set  to  zero,  giving 

Xq  =  [0(  1 , 1(T3)  0(0, 10-3)  0(0,  10“3)  0(1 ,  10~3) 

X0(O,  cxmax10~3)  0(0 ,  czmax10~3)]r  (30) 

where  0(/r ,  cr)  is  a  Gaussian  random  variable  with  mean  /x  and  vari¬ 
ance  <r. 

Finally,  BFGS  can  be  run  for  a  given  number  of  search  directions 
or  until  the  objective  function  value  or  the  magnitude  of  the  gradient 
falls  below  some  tolerance.  Here,  a  set  number  of  search  directions  is 
chosen  with  a  low  number  (15)  for  scaling  and  translation  and  a 
higher  number  (50)  for  affine  transformations  and  point  optimiza¬ 
tion. 

Algorithm  example 

Figure  6  shows  the  progress  of  the  algorithm  through  an  example. 
(Details  of  the  parameters  are  given  in  the  Results  section.)  The 
method  begins  with  a  single  rectangle  (Figure  6a).  The  set  of  points 
that  make  up  the  rectangle  is  labeled  X,  =  {p(1) ,  .  .  .  ,  p(8>}.  Al¬ 
though  the  initial  guess  is  convex,  the  result  (Figure  6a)  is  given  by 

P*  =  CH(X1).  (31) 

The  matrix  of  point  coordinates  P*  is  used  to  compute  the  gravita¬ 
tional  anomaly  (upper  portion  of  Figure  6a)  via  equation  1 1 ,  and  the 
objective  function  value,  given  by  equation  17,  is  0.133.  Initializa¬ 
tion  continues  by  applying  a  few  scale  and  translate  optimizations.  A 
scaling  operation  (equation  6)  is  optimized  first,  given  by 


BFGS  is  applied  to  this  two-parameter  problem,  fg{s„ ,  szz) ,  as  dis¬ 
cussed  in  the  previous  subsection.  As  shown  in  equation  22,  BFGS 
requires  the  gradient  of  the  objective  function.  This  is  obtained  using 
a  finite-difference  approximation;  for  the  example  above,  it  is  given 
by  |- 

fg(.Sxx  "b  h\  ,  Szz)  fg \SXX  ,  ^zz) 


^fg(sxx  >  szz) 


h  t 


fgi'^XX  5  $zz  fgi'^xx  5  '^Z?) 


(33) 


«2 

The  step  size  is  given  by  hk  =  e^chJc“p,  where  emach  is  the  machine 
precision  and  x“p  is  the  expected  magnitude  of  the  kth  argument  as 
discussed  in  the  previous  subsection. 

Note  that  each  gradient  computation  requires  K  +  1  function 
evaluations  if  K  is  the  number  of  arguments.  As  K  grows,  the  finite- 
difference  computations  can  require  many  function  evaluations, 
though,  fortunately,  they  are  all  independent  and  therefore  easily 
parallelizable.  The  3D  version  of  the  algorithm  takes  advantage  of 
this  as  each  function  evaluation  is  more  costly  than  in  the  2D  case. 

Next,  a  translation  optimization  (equation  7)  is  performed,  given 


by 


P*  =  CH(X,)  +R(tx,tz).  (34) 


Again,  this  represents  a  two-parameter  optimization  problem  in  tx 
and  tz.  These  optimizations  are  repeated  three  times,  with  the  final  re¬ 
sult  shown  in  Figure  6b.  The  objective  function  value  at  this  step  is 
0.062. 

Finally,  the  result  is  optimized  using  a  full  affine  transformation 
(equation  5)  to  check  for  convergence.  This  operation  is  given  by 


P*  =  A(jj„,jJCH(X1).  (32) 

(Here  and  below,  the  presentation  is  simplified  by  omitting  the  shift 
of  origin  shown  in  equation  4.) 


P*  =  A (sxx ,  sxz ,  sZJC ,  vJCH( X  0  +  Bfe ,  tz) ,  (35) 

and  the  result  is  shown  in  Figure  6c.  The  objective  function  value  is 
0.04,  which  does  not  meet  the  termination  criteria. 


Figure  6.  An  example  of  the  inversion  process,  (a)  Rectangle,  (b)  Translation  and  scaling  optimizations,  (c)  Affine  optimization,  (d)  Splitting,  (e) 
Affine  optimization,  (f)  Point  optimization  and  final  result. 
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After  the  initialization  stage,  the  algorithm  attempts  to  add  more 
convex  shapes  as  necessary.  The  result  of  the  previous  stage  is  split 
into  two  point  sets,  X  ]  and  X2,  and  optimized  twice,  first  after  split¬ 
ting  at  0°  and  again  at  90°.  Each  separate  split  is  optimized  using  a 
partial  affine  operation  (equation  8),  now  given  by 

P*  =  [AjCs^.i^JCHCXj)  + 

U  [A2(s® ,  4?)CH(X2)  +  B2(ff ,  tf)l  (36) 

which  is  an  eight-parameter  optimization  problem.  In  this  example, 
the  better  choice  was  90°  (Figure  6d),  the  result  of  which  had  an  ob¬ 
jective  function  value  of  0.01.  Again,  the  result  is  optimized  using  a 
full  affine  transformation  to  check  for  convergence,  the  result  of 
which  is  shown  in  Figure  6e.  (The  optimization  problem  here  is  com¬ 
posed  of  12  parameters  and  is  similar  to  that  given  in  equation  36.) 
This  result  (Figure  6e)  has  a  sufficiently  low  error  (1.1  X  10~3)  to 
proceed  to  the  optimize  stage. 

At  this  stage,  two  types  of  optimizations  are  run  until  a  termina¬ 
tion  criterion  is  achieved.  First,  another  full  affine  optimization  is 
performed  (not  shown).  Next,  new  points  can  be  added  to  either  set 
as  described  previously.  In  this  example,  no  additional  points  were 
added  by  the  algorithm.  Finally,  the  point  sets,  X,  =  {pu>, 

.  .  .  ,  p(6,|  and  X2  =  {q(1) ,  .  .  .  ,  q<6)},  are  optimized  directly.  The 
operation  used  by  the  optimizer  is  given  by 

P*  =  CH[X!(p(1),  .  .  .  ,p(6))] 

UCH[X2(q(1\ - q(6))],  (37) 

and  a  total  of  24  parameters  are  required  for  optimization  in  this  case 
(six  points  in  both  X,  and  X2).  The  final  result,  with  an  objective 
function  value  of  3. 8  X  10  s,  is  shown  in  Figure  6f. 

RESULTS 

Several  results,  both  2D  and  3D  and  synthetic  and  measured  are 
presented.  The  first  subsection  gives  results  for  synthetic  magnetic 
anomalies  and  the  second  comprises  synthetic  gravitational  anoma¬ 
lies.  The  final  subsection  gives  results  for  measured  data  sets  of  grav¬ 
itational  data.  Finally,  in  each  2D  figure,  solid  lines  indicate  the  tar- 


Figure  7.  Inversion  of  a  three-mass  structure  using  gravitational 
data.  The  target  is  shown  as  the  solid  line  and  the  inversion  as  the 
dashed  line. 


get  stiucture  and  target  data  (or  previously  published  result)  and 
dashed  lines  indicate  the  results  of  the  inversion. 

Each  example  gives  the  total  number  of  function  evaluations, 
which  is  a  rough  measure  of  the  computation  time.  The  2D  method 
was  implemented  in  MATLAB  and  the  3D  method  was  implement¬ 
ed  in  C+  + .  For  the  MATLAB  version,  typical  run  times  on  a  2-GHz 
Pentium  M  laptop  are  around  7  minutes,  though  it  is  important  to 
note  that  MATLAB  is  interpreted  (not  compiled)  code  and  typically 
is  slower  than  any  compiled  code.  The  3D  version  was  run  on  a  clus¬ 
ter,  with  the  finite-difference  computations  run  in  parallel.  Run  times 
for  this  code  are  on  the  order  of  an  hour. 

Gravitational  anomaly 

As  a  synthetic  example,  a  fault  was  simulated  as  shown  in  the  ex¬ 
ample  given  in  the  previous  section.  Twenty  stations  were  spaced 
equally  over  5  km;  the  density  contrast  was  assumed  to  be  A p 
=  0.276  g/cm3.  Figure  6f  shows  the  target  as  the  solid  line,  the  opti¬ 
mized  result  as  the  dashed  line,  and  the  simulated  and  inverted  gravi¬ 
ty  anomalies  at  the  top.  This  result  required  5486  function  evalua¬ 
tions  to  achieve  an  error  of  3. 8  X  10^5. 

Next,  a  synthetic  structure  made  up  of  three  masses  —  one  large, 
deep  mass  and  two  shallow,  smaller  masses  as  shown  in  Figure  7  — 
was  used  as  a  target.  Again,  a  set  of  25  stations  was  placed  over 
40  km  and  the  density  contrast  was  the  same  as  above.  Figure  7 
shows  the  results  of  the  inversion,  which  achieved  an  error  of  1.9 
X  10~4after  11, 112  function  evaluations. 

Additionally,  the  structure  was  inverted  in  the  presence  of  cor¬ 
rupted  data.  First,  a  noise  level  of  50  dB  was  used,  corresponding  to 
a  noise  level  of  approximately  0.14  mGal  (rms  average).  The  algo¬ 
rithm  was  run  twice,  each  time  with  a  different  set  of  noisy  data.  The 
results  are  shown  in  Figure  8.  Next,  a  noise  level  of  40  dB  was  used, 
giving  an  rms  noise  level  of  approximately  0.48  mGal.  Figure  9 
shows  the  results  for  two  separate  runs. 

As  a  3D  example,  a  two-body  structure  was  simulated  with  two 
randomly  generated  convex  shapes  placed  at  different  depths.  This 
example  used  30  stations  and  profiles  over  a  30  X  30  km  region  and 
a  density  contrast  of  Ap  =  0.5  g/cm3.  Figure  10  shows  the  simulated 
and  inverted  gravity  contours.  The  result  shown  in  Figure  11  re¬ 
quired  12,688  function  evaluations  to  reach  a  final  error  of  1.3 
X  10~3.  In  this  figure  (and  other  3D  plots),  depth  is  indicated  by  the 
shading,  i.e.,  the  shallowest  parts  of  the  structure  are  shown  in  white 
and  the  deepest  are  shown  in  black. 

Magnetic  anomaly 

The  three-mass  structure  from  the  previous  subsection  was  also 
imaged  using  magnetic  data.  The  structure  covers  approximately 
40  km  and  is  approximately  20  km  deep.  A  magnetic  anomaly  was 
generated  assuming  a  susceptibility  of  \  ~  10  ~3  (in  SI  units),  a 
strike  angle  of  60°,  and  an  ambient  magnetic  field  of  50  p,T  at  an  in¬ 
clination  angle  of  0°.  Data  were  measured  at  a  set  of  25  stations  from 
0  to  40  km.  After  27,122  function  evaluations,  the  result  shown  in 
Figure  12  was  achieved;  it  has  an  error  of  2.2  X  10-4. 

Finally,  a  synthetic  3D  example  was  inverted.  The  target  structure 
was  composed  of  three  randomly  generated  convex  shapes  as  shown 
in  Figure  14a.  The  structure  had  a  susceptibility  of  x  =  10'3  and  the 
ambient  field  was  assumed  to  have  a  strength  of  50  p,T  at  an  inclina¬ 
tion  angle  of  60°  and  a  declination  angle  of  0°.  The  final  result,  which 
used  data  from  a  30  X  30  point  grid  over  a  20  X  20  km  area,  had  an 
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Figure  8.  Three-mass  structure  inverted  in  the  presence  of  noisy  data  with  an  S  /N  of  50  dB .  (a)  Run  1 .  (b)  Run  2. 


Figure  9.  Three-mass  structure  inverted  in  the  presence  of  noisy  data  with  an  S  /N  of  40  dB .  (a)  Run  1 .  (b)  Run  2. 


x  (km)  x  (km) 

Figure  1 0.  Gravity  anomaly  contours  (in  milligals)  over  a  synthetic  structure,  (a)  Target,  (b)  Inversion. 
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error  of  1.0  X  10  3  and  required  29,901  function  evaluations.  Re¬ 
sults  of  the  inversion  are  shown  in  Figures  1 3  and  14. 

Measured  data  sets 

Several  measured  data  sets  were  also  inverted,  the  first  of  which  is 
an  aeromagnetic  anomaly  described  in  Murthy  et  al.  (2001).  There,  a 


Figure  1 1 .  (a)  Target  and  (b)  inverted  structures. 


Figure  12.  Inversion  of  a  three-mass  structure  using  magnetic  data. 
Target  is  shown  as  the  solid  line  and  the  inversion  as  the  dashed  line. 


fault  structure  was  assumed  and  several  parameters,  such  as  intensi¬ 
ty  of  magnetization,  depth,  and  position,  were  found.  Here,  an  un¬ 
known,  linearly  varying  regional  magnetic  field  is  assumed  and  the 
two  coefficients  are  found  during  the  optimization.  The  two  coeffi¬ 
cients  of  the  regional  field  were  added  to  the  vector  x  at  each  step. 
The  data  set  comprised  26  stations  over  flat  terrain,  between  0  and 
50  km.  The  values  found  in  Murthy  et  al.  (2001) 
are  used:  Remnant  magnetization  was  assumed 
with  an  intensity  of  magnetization  of  148  nT  at  an 
angle  of  —91°,  a  total  field  inclination  angle  of 
35°,  and  a  strike  angle  of  90°. 

The  final  result,  shown  in  Figure  15,  required 
3403  function  evaluations  to  reach  a  minimum  er¬ 
ror  of  4.2  X  10*3.  The  regional  magnetic  field 
had  a  slope  of  —1.44  nT/km  and  an  offset  of 
66  nT.  Also,  because  a  fault  was  being  imaged, 
the  result  in  Figure  15  actually  extends  to  approx¬ 
imately  300  km,  though  only  the  relevant  portion 
of  the  fault  is  shown.  Good  agreement  is  seen 
with  the  result  in  Murthy  et  al.  (2001),  especially 
the  top  of  the  fault.  Both  results  are  around  7  km. 

The  first  gravity  anomaly,  observed  over  the  Weardale  granite 
body  using  23  stations  across  55  km  of  flat  terrain  and  presented  in 
Murthy  and  Rao  (1993),  is  shown  in  Figure  16.  Again,  a  linearly 
varying  regional  gravitational  field  was  assumed  and  optimized  si¬ 
multaneously  with  the  geometry  at  each  step.  The  optimization  algo¬ 
rithm  required  16,380  function  evaluations  to  achieve  a  final  error  of 
5.0  X  10*3.  Figure  16  shows  the  final  geometry,  with  a  density  con¬ 
trast,  as  assumed  in  Murthy  and  Rao  (1993),  of  Ap  =  —0.13  g/cm3; 
it  is  similar  to  the  result  in  Murthy  and  Rao  (1993).  The  regional  field 
variation  had  a  slope  of  4.72  X  10*3  mgal/km  and  an  offset  of 
—  0. 133  mgal.  As  in  Murthy  and  Rao  (1993),  the  top  of  the  structure 
does  not  outcrop  at  the  surface. 

Next,  a  basin  structure  with  a  parabolic  density  distribution  was 
inverted.  The  Bouguer  gravity  anomaly  data  was  taken  from  Rao  et 
al.  (1994)  and  had  an  initial  density  of  A p0  =  —0.5206  g/cm3  and 
decay  parameter  a  =  0.0576  km* 1 .  The  data  set  was  taken  from  13 
stations  over  48  km  of  flat  terrain.  The  final  result  (Figure  17)  had  an 
error  of  6.0  X  10*4  (5688  function  evaluations)  and  is  compared  to 
the  result  shown  in  Rao  et  al.  (1994)  (solid  line).  As  a  basin,  the  in¬ 
verted  structure  should  have  a  flat  top  across  the  entire  region,  as- 
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Figure  13.  Magnetic  anomaly  contours  (in  nanoteslas)  over  a  structure  composed  of  three  random  bodies,  (a)  Target,  (b)  Inversion. 
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Figure  14.  (a)  Target  and  (b)  inverted  three-body  structure. 


Figure  15.  Inversion  of  aeromagnetic  anomaly  near  Dehri.  Bihar,  In¬ 
dia.  Solution  from  Murthy  et  al.  (2001)  is  shown  as  the  solid  line. 
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Figure  16.  Inversion  of  Weardale  granite  body.  Solution  from  Mur¬ 
thy  and  Rao  (1993)  is  shown  as  the  solid  line. 


Figure  17.  Inversion  of  Los  Angeles  basin.  Solution  from  Rao  et  al. 
(1994)  is  shown  as  the  solid  line. 
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Figure  18.  Gravity  anomaly  contours  (in  milligals)  over  the  Gelibolu  Peninsula,  (a)  Measured,  (b)  Inversion. 


x  (km) 


sumed  in  Rao  et  al.  (1994).  Here,  although  no  assumptions  were 
made,  the  top  of  the  structure  is  mostly  flat  with  some  small  devia¬ 
tion  toward  the  left  side. 

As  a  final  example,  a  measured  3D  data  set  was  inverted.  The  Bou- 
guer  anomaly  was  adapted  from  Figure  12  in  Albora  et  al.  (2007),  as 
shown  in  Figure  18a.  In  this  inversion,  37  stations  were  used  along 
35  profiles  over  a  flat  area  covering  18.5  X  17.5  km.  The  density 
contrast  was  assumed  to  be  A p  =  0.276  g/cm3.  The  final  result, 
whose  anomaly  contour  plot  is  shown  in  Figure  1 8b,  required  42,465 
function  evaluations  and  had  an  error  of  1.56  X  10-2.  Figure  19 
shows  three  views  of  the  interpreted  structure:  top  down  view  (Fig¬ 
ure  19a),  perspective  view  perpendicular  to  the  plane  indicated  by  a 
dashed  line  in  the  previous  figure  (Figure  19b),  and  a  contour  taken 
in  the  plane  indicated  by  the  dashed  line  (Figure  1 9c) .  In  Albora  et  al. 
(2007),  a  fault  structure  was  predicted  along  the  dashed  line  on  the 
right  side  of  Figure  19a.  This  structure  appears  in  the  3D  result 
shown  here,  verified  in  Figure  19c. 


CONCLUSIONS 


x-y  (km) 


Figure  19.  Interpreted  structure  at  the  Gelibolu  Peninsula,  (a)  Top 
view,  (b)  Perspective  view  along  dashed  line,  (c)  Contour  taken  in 
the  plane  indicated  by  the  dashed  line. 


Our  method  for  inverting  gravitational  and  magnetic  anomaly 
data  uses  a  tree-based  geometry  description  that  combines  arbitrary 
convex  shapes  using  Boolean  operations.  This  approach  is  flexible 
in  that  no  knowledge  of  the  number,  approximate  location,  or  com¬ 
plexity  of  the  target  geometry  is  required.  A  BFGS-based  optimiza¬ 
tion  algorithm  uses  successive  runs  of  the  local  optimizer  applied  to 
geometric  transforms  and  actual  points.  Separate  stages  in  the  algo¬ 
rithm  first  determine  approximate  shape  and  complexity  and  finally 
details  of  the  structure.  Results  showed  that  different  types  of  geom¬ 
etries,  both  two  and  three  dimensions,  could  be  reconstructed  accu¬ 
rately  using  both  gravitational  and  magnetic  information. 

In  its  current  form,  the  method  has  some  limitations  that  could  be 
overcome.  First,  only  homogeneous  structures  were  treated  here. 
The  method  could  be  extended  to  inhomogeneous  structures  by  al¬ 
tering  the  operations  of  the  tree  structure.  Material  values  would  be 
added  to  the  leaf  nodes  in  the  tree  and  optimized  along  with  the  ge¬ 
ometry.  Also,  as  geometries  become  more  complex,  the  method  at¬ 
tempts  to  search  larger  spaces.  This  could  be  overcome  by  including 
only  a  few  of  the  total  number  of  leaf  nodes  in  each  optimization.  An¬ 
other  possibility  is  the  use  of  optimizers  other  than  BFGS  to  avoid 
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costly  finite-difference  computations.  Stochastic,  derivative-free 
methods  could  be  used  or  even  automated  differentiation  with  ad¬ 
joint  formulations. 
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